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• ■ Abstract 

O ■ Arithmetic constraints on integer intervals are supported in many con- 

straint programming systems. We study here a number of approaches to 
, implement constraint propagation for these constraints. To describe them 

^ • we introduce integer interval arithmetic. Each approach is explained us- 

\jQ ' ing appropriate proof rules that reduce the variable domains. We compare 

t-H . these approaches using a set of benchmarks. For the most promising ap- 

■ proach we provide results that characterize the effect of constraint prop- 

agation. 

o 

\o . 

1 Introduction 

1.1 Motivation 

> ■ 

• *h . The subject of arithmetic constraints on reals has attracted a great deal of at- 

?\ ' tention in the literature. In contrast, arithmetic constraints on integer intervals 

have not been studied even though they are supported in a number of constraint 
programming systems. In fact, constraint propagation for them is present in 
ECI/PS 6 , SICStus Prolog, GNU Prolog, ILOG Solver and undoubtedly most of 
the systems that support constraint propagation for linear constraints on integer 
intervals. Yet, in contrast to the case of linear constraints — see notably [2] 
— we did not encounter in the literature any analysis of this form of constraint 
propagation. 

In this paper we study these constraints in a systematic way. It turns out 
that in contrast to linear constraints on integer intervals there are a number of 
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natural approaches to constraint propagation for these constraints. They differ 
in the extent to which the constraints are decomposed. 

Even though arithmetic constraints on integer intervals need not be decom- 
posed into atomic arithmetic constraints, as is common practice for constraints 
on reals, we found that it is beneficial to do so: it allows for efficient scheduling 
of the reduction rules and for reuse of auxiliary variables for common subterms 
between constraints. 

It could be argued that since integer arithmetic is a special case of real 
arithmetic, specialized constraint propagation methods for integer arithmetic 
constraints are not needed. Indeed, a constraint satisfaction problem (CSP) 
involving arithmetic constraints on integer variables can be solved using any 
known method for constraints on reals, with additional constraints ensuring 
that the variables assume only integer values. This was suggested in [6] and is 
implemented, for example, in RealPaver [13]. However, a dedicated study and 
implementation of the integer case is beneficial for a number of reasons. 

• In some cases the knowledge that we are dealing with integers yields a 
stronger constraint propagation than the approach through the constraint 
propagation for arithmetic constraints on reals. This can be also bene- 
ficial when we are dealing with hybrid problems that involve arithmetic 
constraints on both integer and real variables. 

• The 'indirect' approach through the reals is based on floating-point num- 
bers, which are of limited precision. This implies that no exact represen- 
tation exists for integers outside certain bounds. We believe that it should 
be possible to deal with large integers precisely, and that we should not 
revert to a floating-point representation when other options exist. Using 
a library like GNU MP [TT] we can use arbitrary length integers (called 
multiple precision integers in GNU MP), whose size is limited only by the 
available memory. 

• Since arithmetic constraints on integer intervals are supported in a number 
of constraint programming systems, it is natural to investigate in a sys- 
tematic way various approaches to their implementation. The approaches 
based on the integers are amenable to a clear theoretical analysis. In par- 
ticular, in Section [9] and Subsection 110. II we provide the characterization 
results that clarify the effect of constraint propagation for the approach 
that emerged in our studies as the fastest. 

An example that supports the first argument is the constraint x ■ y — z, where 
—3 < x < 3, —1 < y < 1, and 1 < z < 2. When all variables are integers, 
there are no solutions having x = 3 or x = —3, and the constraint propagation 
methods that we consider here will actually remove these values from the domain 
of x. However, if these variables are considered to be reals, these values may 
not be removed, and solving the integer problem through constraint propagation 
methods for constraints on reals may lead to a larger search space. 

As an indication that integer representation is not entirely a theoretical 
issue, consider the following benchmark from [6]. Find n integers xi, . . . , x n , 
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CPU time (sec.) 



n 


solutions double 


long long mpz 


13 


22 


0.44 


0.41 1.69 


14 


60 


1.37 


1.35 5.27 


15 


159 


4.45 


4.50 17.44 


16 


377 


14.54 


15.04 57.31 


17 


377 


32.66 


33.54 128.26 


18 


1007 


106.77 


110.98 419.74 



Table 1: Comparison of timing results for various representations 
1 < Xi < n, verifying the conditions 

n n n n 

J\ X i = J\ i > X! < x 2 < ■ ■ ■ < x n . 

i—1 i—1 i—1 i—1 

For n = 10 the initial maximum value of the left-hand side expression of the 
second constraint equals 10 10 , which exceeds 2 32 , the number of values that can 
be represented as 32-bit integers. For n = 16, there is already no signed integer 
representation of this bound in 64 bits. 

To show that arbitrary length integers can be affordable, Table [T] shows 
timing results for three small C++ programs that solve the above benchmark 
via a basic branch-and-propagate search process. These programs differ only in 
the representation of the bounds of the variables, and in the signature of the 
arithmetic operations applied to these bounds: 64-bit floating point numbers 
(double), 64-bit integers (long long), and arbitrary length integers (using the 
mpz data type of the GNU MP library). The programs were compiled using the 
same optimization flags as the default build of the GNU MP library, and the 
reported CPU times are user time in seconds, measured by the time command 
on a 1200 MHz AMD Athlon CPU. 

The results for 64-bit integers and n > 16 could be computed by initializing 
the upper bound of the auxiliary variable equated to the product of all problem 
variables to n!, which works for n < 20. These results indicate that on our hard- 
ware, the 64-bit integer and floating-point implementations are equally efficient, 
while for these specific problem instances, the cost of using arbitrary length in- 
tegers is roughly a factor four. Note that in a full-fledged constraint solver, this 
overhead would be far less prominent, because compared to these small C++ 
programs, a large part of the execution time is spent in the framework that 
coordinates the computation (cf. the results in Subsection 1 11.2[) . 

1.2 Plan of the Paper 

In the next section we provide the relevant background material on CSPs and 
arithmetic constraints. The unifying tool in our analysis is integer interval 
arithmetic that is modeled after the real interval arithmetic (see, e.g., [15]). 
There are, however, essential differences since we deal with integers instead of 



3 



reals. For example, the product of two integer intervals does not need to be 
an integer interval. In Section [3] we introduce integer interval arithmetic and 
establish the basic results. Then in Section 2] we show that using integer interval 
arithmetic we can define succinctly the well-known constraint propagation for 
linear constraints on integer intervals. 

The next three sections, [SJ |5] and \7\ form the main part of the paper. We 
introduce there three approaches to constraint propagation for arithmetic con- 
straints on integer intervals. They differ in the way the constraints are treated: 
either they are left intact, or the multiple occurrences of variables are elimi- 
nated, or the constraints are decomposed into a set of atomic constraints. In 
Section [5] we discuss how these three approaches relate to various methods used 
to solve arithmetic constraints on reals. 

Then in Section [9] we characterize the effect of constraint propagation for 
the atomic constraints. In Section 1101 we discuss in detail our implementation 
of the alternative approaches, and in Section [Tl] we describe the experiments 
that were performed to compare them. They indicate that decomposition of 
the constraints, combined with a scheduling of the reduction rules that respects 
the hierarchical dependencies between the atomic constraints is superior to the 
other approaches. Finally, in Section fTS] we provide the conclusions. 

The preliminary results of this work were reported in [5] and [3] . 

2 Preliminaries 

2.1 Constraint Satisfaction Problems 

We now review the standard concepts of a constraint and of a constraint satis- 
faction problem. Consider a sequence of variables X := x\, . . ., x n where n > 0, 
with respective domains D\, . . ., D n associated with them. So each variable Xi 
ranges over the domain Dj. By a constraint C on X we mean a subset of 
Di x . . . x D n . Given an element d := d%, . . ., dn of D\ x . . . x D n and a sub- 
sequence Y := Xi 1} . . ., Xi l of X we denote by d[Y) the sequence d^ , . . ., di x . In 
particular, for a variable Xi from X, d[xi] denotes di. 

A constraint satisfaction problem, in short CSP, consists of a finite 
sequence of variables X with respective domains T>, together with a finite set C 
of constraints, each on a subsequence of X. We write it as (C ; x% £ D\, . . ., x n € 
D n ), where X := x\, . . ., x n and V := D%, . . ., D n . 

By a solution to (C ; X\ G D\, . . .,x n £ D n ) we mean an element d 6 
Di x . . . x D n such that for each constraint C £ C on a sequence of variables 
X we have d[X] £ C. We call a CSP consistent if it has a solution and 
inconsistent if it does not. Two CSPs with the same sequence of variables 
are called equivalent if they have the same set of solutions. In what follows 
we consider CSPs whose constraints are defined in a simple language and when 
reasoning about them we identify the syntactic description of a constraint with 
its meaning being the set of tuples that satisfy it. 

We view constraint propagation as a process of transforming CSPs that 
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maintains their equivalence. In what follows we define this process by means of 
proof rules that operate on CSPs and preserve equivalence. An interested reader 
can consult pQ or [5] for a precise explanation of this approach to describing 
constraint propagation. 

2.2 Arithmetic Constraints 

To define the arithmetic constraints we use the alphabet that comprises 

• variables, 

• two constants, and 1, 

• the unary minus function symbol ' — ', 

• three binary function symbols, '+','— 'and all written in the infix no- 
tation. 

By an arithmetic expression we mean a term formed in this alphabet and 
by an arithmetic constraint a formula of the form 

s op t, 

where s and t are arithmetic expressions and op € {<,<,=, 7^, >,>}■ For 
example 

x 5 ■ y 2 ■ z 4 + 3x ■ y 3 ■ z 5 < 10 + 4x 4 ■ y 6 ■ z 2 - y 2 ■ x 5 ■ z 4 (1) 

is an arithmetic constraint. Here x 5 is an abbreviation for x ■ x ■ x ■ x ■ x, while 
3x ■ y 3 ■ z 5 is an abbreviation for x ■ y 3 ■ z 5 + x ■ y 3 ■ z 5 + x ■ y 3 ■ z 5 , and similarly 
with the other expressions. If '•' is not used in an arithmetic constraint, we call 
it a linear constraint. 

By an extended arithmetic expression we mean a term formed in the 
above alphabet extended by the unary function symbols and ' ^f-"' for each 
n > 1 and the binary function symbol '/' written in the infix notation. For 

example 

W ■ z 4 )/(a; 2 ■ U 5) (2) 

is an extended arithmetic expression. Here, unlike in |T]), a; 5 is a term obtained 
by applying the function symbol l - 5 ' to the variable x. The extended arithmetic 
expressions will be used only to define constraint propagation for the arithmetic 
constraints. 

Fix now some arbitrary linear ordering -< on the variables of the language. 
By a monomial we mean an integer or a term of the form 



where k > 0, X\, ...,Xk are different variables ordered w.r.t. -<, and a is a non- 
zero integer and rz-i, . . .,nk are positive integers. We call then a;™ 1 • . . . • x^ k the 
power product of this monomial. 
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Next, by a polynomial we mean a term of the form 

where n > 0, at most one monomial rrii is an integer, and the power products 
of the monomials m\, . . .,m n are pairwise different. Finally, by a polynomial 
constraint we mean an arithmetic constraint of the form s op b, where s is a 
polynomial with no monomial being an integer, op G {<, <, =, 7^, >, >}, and b is 
an integer. It is clear that by means of appropriate transformation rules we can 
transform each arithmetic constraint to a polynomial constraint. For example, 
assuming the ordering x -< y -< z on the variables, the arithmetic constraint |T]) 
can be transformed to the polynomial constraint 

2x 5 • y 2 ■ z 4 - 4x 4 • y 6 ■ z 2 + 3x • y 3 ■ z 5 < 10 

So, without loss of generality, from now on we shall limit our attention to 
polynomial constraints. 

Next, let us discuss the domains over which we interpret the arithmetic 
constraints. By an integer interval, or an interval in short, we mean an 
expression of the form 

[a..b] 

where a and b are integers; [a..b] denotes the set of all integers between a and 
b, including a and b. If a > b, we call [a..b] the empty interval and denote it 
by 0. By a range we mean an expression of the form 

x e 1 

where x is a variable and I is an interval. Sets of the form {x € Z\x > a} and 
{x £ Z\x < b} are called extended intervals . 

We link the arithmetic constraints with the notion of a constraint defined 
in the previous section by associating in the standard way with each arith- 
metic constraint its interpretation. For an arithmetic constraint on variables 
Xi, . . .,x n with respective integer interval domains D\,.. .,D n this is a subset 
of Di x ... x D n . 

3 Integer Set Arithmetic 

To reason about the arithmetic constraints we employ a generalization of the 
arithmetic operations to the sets of integers. Here and elsewhere Z, Af, and TZ 
denote the sets of all integers, natural numbers, and reals, respectively. 

3.1 Definitions 

For X, Y sets of integers we define the following operations: 
• addition: 

X + Y := {x + y | x G X,y £ Y}, 
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• subtraction: 

X-Y:={x-y\xeX,yeY}, 

• multiplication: 

X-Y :={x-y\xeX,yeY}, 

• division: 

X/Y :={ueZ\3xe X3y eYu-y = x}, 

• exponentiation: 

X n := {x n \ xeX}, 
for each natural number n > 0, 

• root extraction: 

Va := {x e z | x n e X}, 

for each natural number n > 0. 

All the operations except division and root extraction are defined in the 
expected way. We shall return to the division operation in Section At the 
moment it suffices to note the division operation is defined for all sets of integers, 
including Y — and Y = {0}. This division operation corresponds to the 
following division operation on the sets of reals introduced in [T5] : 

X Y := {u e Tl | 3x e X3y eYu-y = x}. 

For an integer or real number a and op G {+, — , •, /, 0} we identify a op X with 
{a} op X and X op a with X op {a}. 

To present the rules we are interested in we shall also use the addition and 
division operations on the sets of reals. Addition is defined in the same way as 
for the sets of integers, and for division we use the operator defined above. In 
|15j it is explained how to implement these operations on, possibly unbounded, 
real intervals. 

Further, given a set A of integers or reals, we define 

- A := {x e Z \ 3a £ A x < a}, 

- A := {x E Z | 3a e A x > a}, 

so for example -M = Z, and -{— 1, 1} and -(—2, 2) both denote the extended 
interval of all integers greater than or equal to —1, where (—2,2) denotes an 
open interval of real numbers. 

When limiting our attention to intervals of integers the following simple 
observation is of importance. 

Note 3.1 For X,Y integer intervals and a an integer the following holds: 

• XC\Y,X + Y, X — Y are integer intervals. 
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• X/{a} is an integer interval. 

• X Y does not have to be an integer interval, even if X = {a} orY = {a}. 

• X/Y does not have to be an integer interval. 

• For each n > 1, X n does not have to be an integer interval. 

• For odd n > 1, \/~X is an integer interval. 

• For even n > 1, \/~X is an integer interval or a disjoint union of two 
integer intervals. □ 

For example in the following cases we get intervals as outcomes: 

[2..4] + [3..8] = [5..12], 

[3..7] - [1..8] = [-S..6], 



S/[-30..100] = [-3..4], 



{/[-100..9] = [-3..3], 
while in the following ones not: 

[3..3]-[1..2]={3,6}, 

[3..5]/[-1..2] = {-5,-4,-3,2,3,4,5}, 
[-3..5]/[-1..2] = Z, 
[1..2] 2 = {1,4}, 



^[119] = [-3.. - 1] U [1..3]. 

To deal with the problem that non-interval domains can be produced by some 
of the operations we introduce the following operation on the sets of integers: 

. J smallest integer interval containing X if X is finite, 

yZ otherwise. 

For example int([3..5]/[-1..2]) = [-5. .5] and int([-3..5]/[-1..2]) = Z. 



3.2 Implementation 

To define constraint propagation for the arithmetic constraints on integer in- 
tervals we shall use the integer set arithmetic, mainly limited to the integer 
intervals. This brings us to the discussion of how to implement the introduced 
operations on the integer intervals. Since we are only interested in maintaining 
the property that the sets remain integer intervals or the set of integers Z we 
shall clarify how to implement the intersection, addition, subtraction and root 
extraction operations of the integer intervals and the int(-) closure of the multi- 
plication, division and exponentiation operations on the integer intervals. The 
case when one of the intervals is empty is easy to deal with. So we assume that 
we deal with non-empty intervals [a..b] and [c.d], that is a < b and c < d. 
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Intersection, addition and subtraction It is easy to see that 



[a..b] n [c.d] = [max(a, c)..min(6, d)], 

[a..b] + [c.d] = [a + c .. b + d], 
[a..b] — [c.d] = [a - 



c .. b - 
d..b- 



So the interval intersection, addition, and subtraction are straightforward to 
implement. 

Root extraction The outcome of the root extraction operator applied to 
an integer interval will be an integer interval or a disjoint union of two integer 
intervals. We shall explain in Section[5]why it is advantageous not to apply int(-) 
to the outcome. This operator can be implemented by means of the following 
case analysis. 



Case 1. Suppose n is odd. Then 



Case 2. Suppose n is even and b < 0. Then 



Vb 



Case 3. Suppose n is even and b > 0. Then 

v / Ml=[- [l ^|J..-D v^Tlllu 

where a + := max(0, a). 



Multiplication For the remaining operations we only need to explain how to 
implement the int(-) closure of the outcome. First note that 

int([o..&] • [c.d]) = [min(A)..max(A)], 

where A — {a ■ c, a ■ d, b ■ c, b • d}. 

Using an appropriate case analysis we can actually compute the bounds of 
int([a..fe] • [c.d]) directly in terms of the bounds of the constituent intervals. 



Division In contrast, the int(-) closure of the interval division is not so straight- 
forward to compute. The reason is that, as we shall see in a moment, we cannot 
express the result in terms of some simple operations on the interval bounds. 

Consider non-empty integer intervals [a..b] and [c.d]. In analyzing the out- 
come of int([a..6]/[c..d]) we distinguish the following cases. 

Case 1. Suppose € [a..b] and G [c.d]. 

Then by definition int([a..6]/[c..d]) = Z. For example, 

int([-1..100]/[-2..8]) = Z. 
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Case 2. Suppose ^ [a..b] and c = d = 0. 

Then by definition int([a..6]/[c..ef]) = 0. For example, 

int([10..100]/[0..0]) = 0. 

Case 3. Suppose ^ [a..b] and c < and < d. 
It is easy to see that then 

int([a..6]/[c..d]) = [-e..e], 

where e = max(|a|, \b\). For example, 

int([-100.. - 10]/[-2..5]) = [-100. .100]. 

Case 4- Suppose ^ [a..b] and either c = and d ^ or c ^ and d = 0. 
Then int([a..6]/[c..d]) = int([a..b]/([c..d] - {0})). For example 

int([1..100]/[-7..0]) = int([1..100]/[-7.. - 1]). 

This allows us to reduce this case to Case 5 below. 
Case 5. Suppose ^ [c.d]. 

This is the only case when we need to compute int([a..6]/[c..d]) indirectly. 
First, observe that we have 

mt([a..b]/[c..d]) C [\mm(A)] .. Lmax(A)J], 

where A — {a/c, a/d, b/c, b/d}. 

However, the equality does not need to hold here. Indeed, note for example 
that int([155..161]/[9..11]) = [16. .16], whereas for A = {155/9,155/11,161/9, 
161/11} we have [min(y4)] = 15 and [max(A)J — 17. The problem is that the 
value 16 is obtained by dividing 160 by 10 and none of these two values is an 
interval bound. 

This complication can be solved by preprocessing the interval [c.d] so that 
its bounds are actual divisors of an element of [a..b]. First, we look for the least 
c! E [c.d] such that 3x £ [a..b] 3u £ Z u-c' — x. Using a case analysis, the latter 
property can be established without search. Suppose for example that a > and 
c > 0. In this case, if c' • [b/c'\ > a, then c' has the required property. Similarly, 
we look for the largest d! £ [c.d] for which an analogous condition holds. Now 
mt([a..b]/[c.d]) = [\min(A)]..[max(A)\], where A = {a/d, a/d', b/c', b/d'}. 

In view of this auxiliary computation (in case when ^ [c.d]) we shall 
introduce in Section [10] a modified division operation with a more direct imple- 
mentation. 

Exponentiation The int(-) closure of the interval exponentiation is straight- 
forward to implement by distinguishing the following cases. 

Case 1. Suppose n is odd. Then 

int([a..&n = [a n ..b n ]. 
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Case 2. Suppose n is even and < a. Then 

int([a..&n = [a n ..b n ]. 
Case 3. Suppose n is even and b < 0. Then 

int([a..&n = [b n ..a n ]. 
Case 4- Suppose n is even and a < and < b. Then 
int([a..6] n ) = [0.. max(a", &")]. 

3.3 Correctness Lemma 

Given now an extended arithmetic expression s each variable of which ranges 
over an integer interval, we define int(s) as the integer interval or the set Z 
obtained by systematically replacing each function symbol by the application 
of the int(-) operation to the corresponding integer set operation. For example, 
for the extended arithmetic expression s := %J (y 2 ■ z 4 )/(x 2 ■ u 5 ) of |2]) we have 

int(s) = int(^/int(int(y 2 ) • int(Z 4 ))/ int(int(X 2 ) • int(?7 5 ))), 

where we assume that x ranges over X, etc. 

The discussion in the previous subsection shows how to compute int(s) given 
an extended arithmetic expression s and the integer interval domains of its 
variables. 

The following lemma is crucial for our considerations. It is a counterpart 
of the so-called 'Fundamental Theorem of Interval Arithmetic' established in 
|17j . Because we deal here with the integer domains an additional assumption 
is needed to establish the desired conclusion. 

Lemma 3.2 (Correctness) Let s be an extended arithmetic expression with 
the variables X\, . . .,x n . Assume that each variable Xi of s ranges over an integer 
interval Xj. Choose S Xi for i G [l..n] and denote by s(a±, . . ., a n ) the result 
of replacing in s each occurrence of a variable Xi by a.j . 

Suppose that each subexpression of s(a\, . . ., a n ) evaluates to an integer. 
Then the result of evaluating s(a\, . . ., a n ) is an element o/int(s). 

Proof. The proof follows by a straightforward induction on the structure of s. 

□ 

4 An Intermezzo: Constraint Propagation for 
Linear Constraints 

Even though we focus here on arithmetic constraints on integer intervals, it 
is helpful to realize that the integer interval arithmetic is also useful to define 
in a succinct way the well-known rules for constraint propagation for linear 



11 



constraints (studied in detail in [II]). To this end consider first a constraint 
S" =1 aj • Xi = b, where n > 0, 01, . . .,a n are non-zero integers, Xi, . . .,x n are 
different variables, and b is an integer. We rewrite this constraint n times, 
each time isolating one variable, to obtain an extended arithmetic expression 
for each variable Xj . To each of these extended arithmetic expressions we apply 
then the int operation of Subsection [231 which yields an update for the domain 
of the corresponding variable Xj. To reason about this procedure we can use 
the following rule parametrized by j G 

LINEAR EQUALITY 

(gjLjfflj ' x i = b j xi € . . ., Xj G Dj, . . ., x n G D„) 
(S" =1 aj ■ Xi = b ; xi G Di, . . .,Xj G Dj, . . .,x„ G D n ) 

where Dj := Dj n int ((& - S ie [i..„]-{ : ,}a l • Xi)/a,J. 
Note that by virtue of Note 13.11 

D j = D o n ( & ~ s ie[i..n]-{j} int K ' D i))/aj- 

To see that this rule preserves equivalence, first note that all our reduction 
rules compute the domain updates via intersection with the original domain, 
preventing that domains are extended by their application. Further, suppose 
that for some d\ G D\, . , ., d n G D n we have £™ =1 ai • di = b. Then for j G [l..n] 
we have 

^' = Q 3 ~ ^i£[l..n]-{j}«i ' di)/aj 

which by the Correctness Lemma 13.21 implies that 

dj G int nb- £,' e [i.. n ]_{j}ffli -x^/aj), 

i.e., dj G Dj. 

Next, consider a constraint E" =1 a^ • Xj < 6, where oi, . . ., a ra ,xi, . . .,x„ and 
6 are as above. To reason about it we can use the following rule parametrized 
by j £ [l..n]: 

LINEAR INEQUALITY 

(gjLjfflj ' Xj < b ; xi € Di, . . .,Xj € Dj, . . .,x n £ At) 
(S" =1 ai ■ Xi < b ; xi G D\, . . .,Xj G Dj, . . .,x„ G D n ) 

where Dj := Dj n (-int (6 - £ ie [i.. n ]_{j}a, • x,-)/aj). 

To see that this rule preserves equivalence, suppose that for some d\ G 
Di, . . ., g?„ G D n we have S^_ 1 a, ■ di < b. Then <Zj • < 6 — S ie [ 1 _ n ]_y-j.aj • cZj. 
By the Correctness Lemma I3~^1 

b — ^ie[i..n]-{j}0-i ■ di G int(6 — S jg [x..n]— y}^* ' 
so by definition 

dj • <ij G- int (6 - E, e [i..„]_{j}ai • x») 
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and consequently 

dj £- int(6 - ^i & [x.. n ]-{j}ai ■ Xi)/a r 

This implies that dj £ DL 

As an alternative to evaluating int(Sjg[i .. n i_/j\aj on every application of 
the LINEAR EQUALITY and LINEAR INEQUALITY rules, we could main- 
tain the interval int(S™ai • x{) in an auxiliary variable, and subtract int(aj • Xj) 
from it. This corresponds to the two-step propagation described in [14]. If 
changes to Dj are propagated back to the auxiliary variable, this does not af- 
fect the reduction achieved by the subsequent applications of the rules, while 
the number of interval arithmetic operations involved in the application of a 
rule becomes constant, instead of linear in the number n of variables. 



5 Constraint Propagation: Direct Approach 

As a first approach to constraint propagation for arithmetic constraints on inte- 
ger intervals, we propose to use the constraints directly, in their original form. 
This is an extension of the approach of Section [4] from linear constraints to 
general arithmetic constraints, and entails that these constraints are rewritten 
to isolate all occurrences of each variable. The resulting extended arithmetic 
expressions are then evaluated to obtain updates for the isolated variables. 
The following example illustrates this approach. Consider the constraint 

x 3 ■ y - x < 40 

and the ranges x £ [1..100] and y £ [1..100]. We can rewrite it as 

x < 



^(40 + x)/y\ (3) 

since x assumes integer values. The maximum value the expression on the 
right-hand side can take is [v 7 !!)], so we conclude x < 5. By reusing ([3]), now 
with the information that x £ [1..5], we conclude that the maximum value the 
expression on the right-hand side of ((3|) can take is actually [s/45] , from which 
it follows that x < 3. 

In the case of y we can isolate it by rewriting the original constraint as 
y < 40/a; 3 + l/x 2 from which it follows that y < 41, since by assumption x > 1. 
So we could reduce the domain of x to [1..3] and the domain of y to [1..41]. This 
interval reduction is optimal, since x = l,y — 41 and x = 3, y = 1 are both 
solutions to the original constraint x 3 ■ y — x < 40. So rewriting the constraint 
as x > x 3 ■ y — 40 does not yield a new lower bound for x. 

More formally, consider a polynomial constraint E™ x m,i = b where m > 0, 
no monomial m, is an integer, the power products of the monomials are pairwise 
different, and b is an integer. Suppose that x\, . . ., x n are its variables ordered 
w.r.t. -n. 
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Select a non-integer monomial mi and assume it is of the form 



«■!!: ■■■■■Hi- 

where k > 0, y\,...,yk are different variables ordered w.r.t. -<, a is a non- 
zero integer and m,...,nk are positive integers. So each yi variable equals to 
some variable in {xi, . . .,x n }. Suppose that y p equals to Xj. We introduce the 
following proof rule: 

POLYNOMIAL EQUALITY 
(E" =1 mi = b ; x\ G D u . . .,Xj e Dj, . . .,x n e D n ) 



(E™ =1 m 4 = b ; x\ G D\, . . .,Xj G D'-, . . ., x n G D. 



where Dj := int n y int ((6 — S, e [i .. m ]_{/}rrii)/s) J , with s := a - y" 1 • . . . • 

Vp-i ' Vp+i ■••'Vk ■ 

To see that this rule preserves equivalence, choose some d\ £ D%, . . ., dn G 
D n . To simplify the notation, given an extended arithmetic expression t denote 
by t' the result of evaluating t after each occurrence of a variable Xi is replaced 
by di. 

Suppose that E™ ^ = 6. Then 

d™ p ■ s' = b- £ 4e [i.. m ]-{!} m i: 
so by the Correctness Lemma 13.21 applied to b — £j 6 h, ro ]_mmj and to s 
d r - p e mt(b- E ie [ 1 .. m ]_ W Jn i )/int(s). 

Hence 



dj G "y/hrtO - E ie[ i.. m ]_{ i }m i )/ int(s) 
and consequently 

dj G int \ Dj n "^/int ((6 - E ie[1 .. m] _ {i} mi)/s) ^ 



i.e., G D'j 



Note that we do not apply int(-) to the outcome of the root extraction 
operation. For even n p this means that the second operand of the intersection 
can be a union of two intervals, instead of a single interval. To see why this is 
desirable, consider the constraint x 2 — y = in the presence of ranges x € [0..10], 
y G [25. .100]. Using the int(-) closure of the root extraction we would not be 
able to update the lower bound of a; to 5. 

Next, consider a polynomial constraint E™ 1 rrii < b. Below we adopt the 
assumptions and notation used when defining the POLYNOMIAL EQUALITY 
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rule. To formulate the appropriate rule we stipulate that for the extended 
arithmetic expressions s and t 

int((^«)/t) := -Qn -Q, 

with Q = (^int(s))/int(t). 

To reason about this constraint we use the following rule: 

POLYNOMIAL INEQUALITY 

(Sg^mj < b ; x x g £>i,... ,Zj G Dj,...,x n G £>„) 
(E" =1 mi < 6 ; xi € £>i, . . .,Xj G . ..,x n £ D n ) 

where £K := hit n "0nt (^(6 - £i G [i.. m ]_{;}mi)/s) ) , with s := a ■ y" 1 • 

■ ■ ■ • 2/p_i • Vp+i ■ ■ ■ ' Vk ■ 

To prove that this rule preserves equivalence, choose some d\ G Dx,...,d n G 
D n . As above given an extended arithmetic expression t we denote by t 1 the 
result of evaluating t when each occurrence of a variable Xi in t is replaced by 
di. 

Suppose that S^m^ < 6. Then 

d™ p ■ s' <b - E ie [i.. m ]_{i}m-. 
By the Correctness Lemma [3T21 

& _ ^ie[l..m]-{i} m i e mt (^ _ ^i£[l..m]-{l} m i), 

so by definition 

d" p • s' G- int(6 - S ie[1 .. m] _ w mi). 
Hence by the definition of the division operation on the sets of integers 

d'- p S- int(&- E. ie[1 .. m] _ {z} m J )/int(s) 

Consequently 

G "^/- int(6 - E ie [i.. m ]-{z}m i )/ int(s). 

This implies that G D'-. 

Note that the set - int(6 — £j e n m i_.nimi), which occurs when the expres- 
sion for Dj is expanded according to the above definition of int((— s)/t), is not 
an interval. So to properly implement this rule we need to extend the imple- 
mentation of the division operation discussed in Subsection l3.2l to the case when 
the numerator is an extended interval. 

If the sum of the intervals associated with each of the monomials in a poly- 
nomial constraint is maintained in an auxiliary variable, as we discussed at the 
end of Section[4]for linear constraints, then the rules can be applied using a con- 
stant number of interval additions. However, interval division is not the inverse 
operation of interval multiplication, so the same technique cannot be applied 
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to monomials, and the number of multiplications and exponentiations will be 
linear in the size of the monomial s. 

In an optimized version of the direct approach we simplify the fractions of 
two polynomials by splitting the division over addition and subtraction and by 
dividing out common powers of variables and greatest common divisors of the 
constant factors. Subsequently, fractions whose denominators have identical 
power products are added. We used this optimization in the initial example by 
simplifying (40 + x)/x 3 to 40/x 3 + 1/x 2 . The reader may check that without 
this simplification step we can only deduce that y < 43. 

To provide details of this optimization, given two monomials s and t, we 
denote by 

[s/t] 

the result of performing this simplification operation on s and t. For example, 
[(2 -a; 3 ■ y) / '(4 ■ x 2 )] equals (x-y)/2, whereas [(4 • x 3 -y)/(2-y 2 )] equals (2-x 3 )/y. 

Because the validity of the simplification depends on the sign of the denom- 
inator, we assume that the domains of the variables y\, . . ., y p — i, y p +i, . . ., y n of 
mi do not contain 0. For a monomial s involving variables ranging over the in- 
teger intervals that do not contain 0, the set int(s) either contains only positive 
numbers or only negative numbers. In the first case we write sign(s) = + and 
in the second case we write sign(s) = — . 

The new domain of the variable xj in the POLYNOMIAL INEQUALITY 
rule is defined using two sequences m' Q ...m' n and s' ...s' n of extended arithmetic 
expressions such that 

m' /s' Q — [b/s] and rn! i j s\ = —[rrii/s] for i £ [l..m]. 

Let S := {s^ | i € [0..m] — {I}} and for an extended arithmetic expression f G S 
let I t :— {i E [Q..m] — {1} | — t}. We denote then by p t the polynomial 
Sie/ m 'i- The new domains are then defined by 



D'j := int [Dj n < int (S teS Pt t) J 

if sign(s) = +, and by 

D' 3 := int \Dj n Y> int(E teS p t t) J 

if sign(s) = — . Here the int(s) notation used in the Correctness Lemma \3. 21 is 
extended to expressions involving the division operator on real intervals in 
the obvious way. We define the int(-) operator applied to a bounded set of reals, 
as produced by the division and addition operators in the above two expressions 
for Dp to denote the smallest interval of reals containing that set. 

Returning again to the discussion of the two-step propagation technique 
of [2] , which we started at the end of Section [H note that in this case, the 
int(-) operation is applied after removing the common power products. For this 
reason, there is no straightforward way to calculate int(£ te s pt t) from the 
sum of all intervals associated with the monomials of a polynomial constraint. 
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6 Constraint Propagation: Partial Decomposi- 
tion 



As a second approach to constraint propagation for arithmetic constraints on 
integer intervals, we limit our attention to a special type of polynomial con- 
straints, namely the ones of the form s op 6, where s is a polynomial in which 
each variable occurs at most once and where b is an integer. We call such a con- 
straint a simple polynomial constraint. By introducing auxiliary variables 
that are equated with appropriate monomials we can decompose any polyno- 
mial constraint into a sequence of simple polynomial constraints. This allows us 
also to compute the integer interval domains of the auxiliary variable from the 
integer interval domains of the original variables. We apply then to the simple 
polynomial constraints the rules introduced in the previous section. 

To see that the restriction to simple polynomial constraints can make a 
difference consider the constraint 

100a; • y - lOy ■ z = 212, 

and ranges x, y, z e [1..9]. We rewrite it into the sequence 

u — x ■ y, v — y ■ z, lOOu — 10i> = 212, 

where u,v are auxiliary variables, each with the domain [1..81]. 

It is easy to check that the POLYNOMIAL EQUALITY rule introduced in 
the previous section does not yield any domain reduction when applied to the 
original constraint 100a; • y — lOy ■ z — 212. In the presence of the discussed 
optimization the domain of x gets reduced to [1..3]. 

However, if we repeatedly apply the POLYNOMIAL EQUALITY rule to 
the simple polynomial constraint lOOtt — lOv — 212, we eventually reduce the 
domain of u to the empty set (since this constraint has no integer solution 
in the ranges u, v S [1..81]) and consequently can conclude that the original 
constraint 100a; • y — 10y • z — 212 has no solution in the ranges x,y, z € [1..9], 
without performing any search. Note that this effect still occurs if we replace 
one occurrence of y by a fresh variable with the same domain. 

As noted in [S], decomposing constraints also prevents the evaluation of 
subexpressions whose domains did not change, which may reduce the number 
of interval arithmetic operations performed during constraint propagation. In 
our case duplicate occurrences of variables are removed, so the reduction rules 
additionally become idempotent. However, this can be seen as a side-effect: 
rules still update variables that they depend on, only now this update is indirect, 
through other variables. 

Consider for example the constraint a; 3 • y — x < 40 of Section [5l If we 
rewrite this constraint as u — x < 40, with u — x 3 ■ y and x,y £ [1..100], then 
via u < x + 40 we can set the upper bound for u to 140. Via x = yujy we can 
then set the upper for a; to 5. This allows us to set the upper bound for u to 
45 via u < x + 40, etc. From this point of view the auxiliary variables, and the 
idempotence that they entail, can be seen as an optimization that prevents the 
evaluation of expressions that will not lead to further domain updates. 
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7 Constraint Propagation: Full Decomposition 



In this third approach we focus on a small set of 'atomic' arithmetic constraints. 
We call an arithmetic constraint atomic if it is in one of the following two forms: 

• a linear constraint, 

• x ■ y = z. 

Using appropriate transformation rules involving auxiliary variables we can 
decompose any arithmetic constraint to a sequence of atomic arithmetic con- 
straints, similar to the decomposition of linear constraints into constraints on 
groups of three variables in clp(FD) S\. In this transformation, as with partial 
decomposition, the auxiliary variables are equated with monomials, so we can 
easily compute their domains. 

We explained already in Section[4]how to reason about linear constraints. For 
a treatment of disequalities see, e.g., [TH I2U] . Next, we focus on the reasoning 
for the multiplication constraint x ■ y — z in the presence of the non-empty 
ranges x £ D x , y £ D y and z £ D z . To this end we introduce the following 
three domain reduction rules: 

MULTIPLICATION 1 



(x 


y 


= z ; x £ D x ,y £ D y , z 


G D z ) 


{x 


I) 


— z ; x £ D x , y £ D y , z 


G^) 






MULTIPLICATION 2 




(x 


y 


= z ; x £ D x ,y £ D y ,z 


G D z ) 


(x 


V 


= z ; x £ D' x ,y £ D y ,z 


G D z ) 






MULTIPLICATION 3 




(x 


y 


= z ; x £ D x ,y £ D y ,z 


eD z ) 


(x 


V 


= z ; x £ D x ,y £ D' y ,z 


G D z ) 



where D' z := D z n int(L> x • D y ), D' x := D x n \nt{D z /D y ), and D' y := D y n 
mt(D z /D x ). 

The way we defined the multiplication and the division of the integer inter- 
vals ensures that the MULTIPLICATION rules 1,2, and 3 are equivalence pre- 
serving. Consider for example the MULTIPLICATION 2 rule. Take some a £ 
D X1 b £ D y and c £ D z such that a-b = c. Then a £ {x £ Z \ 3z £ D z 3y £ D y x ■ y 
so a £ D z /D y and a fortiori a £ u&(D z / D y ). Consequently a £ D x flint (D z /D y ). 
Because we also have D x n mt(D z /D y ) C D x , this shows that the MULTIPLI- 
CATION 2 rule is equivalence preserving. 

The following example shows an interaction between all three MULTIPLI- 
CATION rules. 
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Example 7.1 Consider the CSP 

{x-y = z; x£ [1-20], ye [9..11],ze [155.. 161]). (4) 

To facilitate the reading we underline the modified domains. An application 
of the MULTIPLICATION 2 rule yields 

(x-y = z; x£ [16.. 16] , y £ [9. .11], z £ [155. .161]) 

since, as already noted in Subsection 13. 2[ [155..161]/[9..11]) = [16. .16], and 
[1..20] n int([16..16]) = [16.. 16]. Applying now the MULTIPLICATION 3 rule 
we obtain 

(x-y = z; x£ [16..16],y £ [10-10] , z £ [155. .161]) 

since [155..161]/[16..16] = [10.. 10] and [9..11] n int([10..10]) = [10.. 10]. Next, by 
the application of the MULTIPLICATION 1 rule we obtain 

(x-y = z ; x£ [16-16], y£ [10-10], z£ [160-160]) 

since [16. .16] • [10.. 10] = [160.. 160] and [155. .161] H int([160..160]) = [160.. 160]. 
So using all three multiplication rules we could solve the CSP (UJ). □ 

Now let us clarify why we did not define the division of the sets of integers 
Z and Y by 

Z/Y :={z/y£Z\y£Y,z£Z,y^0}. 

The reason is that in that case for any set of integers Z we would have Z/{0} = 
0. Consequently, if we adopted this definition of the division of the integer 
intervals, the resulting MULTIPLICATION 2 and 3 rules would not be anymore 
equivalence preserving. Indeed, consider the CSP 

(x-y = z; x£ [-2-1], y £ [0-0], z £ [-8.. 10]). 

Then we would have [-8..10]/[0..0] = and consequently by the MULTIPLI- 
CATION 2 rule we could conclude 

(x-y = z; x£(b,y£ [0-0], z£ [-8. .10]). 

So we reached an inconsistent CSP while the original CSP is consistent. 

The transformation to atomic constraints can strengthen the reduction. Con- 
sider for example the simple constraint 

w • x • y ■ z = 24 

with w = 4 and x, y, z £ [1-4]. Application of the POLYNOMIAL EQUALITY 
rule does not reduce any of the domains, but if we replace the constraint with 

u ■ v — t, w ■ x = u, y ■ z = v 
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with additional ranges t — 24, u £ [4. .16], and v £ [1..16], then by application 
of the MULTIPLICATION 3 rule to u ■ v = t we reduce the domain of v to 
[2. .6]. Next, by application of the MULTIPLICATION 2 rule to the same 
constraint we reduce the domain of u to [4.. 12], and finally by application of the 
MULTIPLICATION 3 rule to w ■ x = u we reduce the domain of x from [1..4] 
to [1..3]. Note, however, that this effect depends on the decomposition. If we 
had decomposed the constraint as 

z ■ (y ■ ( x ■ w)) = 24 

with an auxiliary variable introduced for each pair of matching brackets, then 
we would not have been able to reduce any of the domains of x, y, and z. 

In the remainder of the paper we will also consider variants of the full de- 
composition approach where we allow squaring and exponentiation as atomic 
constraints. For this purpose we explain the reasoning for the constraint x = y n 
in the presence of the non-empty ranges x € D x and y £ D y , and for n > 1. To 
this end we introduce the following two rules: 

EXPONENTIA TION 

(xj=y™j x € D x ,y € D y ) 
(x = y n ; x£D' x ,y£D y ) 

ROOT EXTRACTION 

(x = y n ; x£D x ,y£ D y ) 
(x = y n ; x£D x ,yeD' y ) 

where D' x := D x n int (£>£), and D' y := iat(D y H \[LT X ). 

To prove that these rules are equivalence preserving suppose that for some 
a 6 D x and b 6 D y we have a — b n . Then a £ D y , so a £ int(D™) and conse- 
quently a £ D x n 'mt(Dy). Also b £ \/D x , so b £ D y n \/D x , and consequently 
b £ hit (A, n y/L%). 

With exponentiation as an atomic constraint, full decomposition leads to 
idempotent rules, and the discussion at the end of Section [6] applies. 



8 Relation to Hull and Box Consistency 

In this section we relate the three approaches introduced above to the well- 
known methods for constraint propagation of arithmetic constraints on real vari- 
ables, whose domains are represented by floating-point intervals. An overview 
of these methods is provided in [9] . Floating-point intervals are intervals of 
reals, with bounds from a finite set T C TZ U {— oo, oo} of floating-point num- 
bers that contains representations — oo and oo for plus and minus infinity. For 
floating-point intervals, the counterpart of the int(-) operation is the hull of 
a set of real numbers defined as the smallest floating-point interval containing 
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the set. Ideally, for an arithmetic constraint c on the variables xi, . . ., x n with 
respective floating-point interval domains D%, . . .,D n we would like to enforce 
hull consistency , which entails that for all i S [l..n] 

Di = hull(xj £ TZ | 3xi G Di, . . . ,Xi-i € A-i,x t+ i € A+i, ■ ■ ■ , x„ € D„ 

(xi, . . . ,x n ) e c). 

However, no efficient procedure exists for enforcing hull consistency on ar- 
bitrary arithmetic constraints. Therefore, the natural approach is to first de- 
compose constraints into atomic constraints, each containing a single arithmetic 
operation. Maintaining hull consistency for the decomposed constraints can be 
done efficiently, using proof rules similar to the ones that we introduced, but hull 
consistency for the resulting decomposed CSP is a weaker notion of consistency 
than hull consistency for the original CSP. 

Our full decomposition approach can be seen as the integer interval equiva- 
lent of the method for computing hull consistency for a decomposition that we 
just described, with the exception that linear constraints are not decomposed 
further. In the floating-point case, because of the accumulation of the round- 
ing errors, the characterization of the resulting form of constraint propagation 
in terms of hull consistency is possible only if all constraints, including linear 
constraints, are decomposed into single-operator constraints. 

To illustrate this consider the constraint x + y + z — w with the variables 
ranging over the floating-point intervals D x , D y , D z and D w . When we evalu- 
ate D x + D y + D z using the floating-point interval arithmetic to compute an 
update for D w , we have three options which two intervals to add first. Because 
the floating-point addition is non-associative, we actually compute the hull of a 
decomposition that has an extra variable added for either x + y, x + z or y + z, 
and the resulting interval is potentially a proper superset of \m\\(D x + D y + D z ). 
Moreover, different rewritings of the constraint correspond to different decom- 
positions, and although this need not be a problem in practice, the resulting 
form of local consistency is no longer clearly defined. 

In contrast, for integer intervals, we do not need to deal with the accumu- 
lation of the rounding errors and the linear constraints can be left intact. Our 
other two approaches can be seen as variants of the full decomposition approach 
that exploit this property further: for partial decomposition we allow more than 
one multiplication per proof rule, and in the direct approach the decomposition 
is not made explicit at all. Apart from these variations, all three approaches 
are the same in one important aspect: multiple occurrences of the same vari- 
able are treated as different variables. To illustrate this, consider the constraint 
x 3 + x = 0, with x € [— 1..1]. While x = is the unique solution, none of our 
three approaches will be able to reduce the domain of x. The reason is that 
the two occurrences of x are essentially treated as different variables in the re- 
duction rules. This problem is known as the dependency problem of interval 
arithmetic. 

In the context of constraints on reals [5] proposed to deal with the depen- 
dency problem using the notion of box consistency. It is a weaker notion of local 
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consistency than hull consistency, but is potentially stronger than hull consis- 
tency for the decomposition of a constraint into atomic constraints (see, e.g., 
[9 ). Enforcing box consistency, as described in [3T], consists of fixing the do- 
mains of all variables except one, and then narrowing the domain of this variable 
by iteratively instantiating it with subdomains at the boundary of the original 
domain, each time verifying consistency of the constraint in the presence of the 
domains of the other variables, and subtracting the subdomain from the original 
domain if the instantiation leads to a failure. 

The second step of the 'trial-and-prune' procedure for enforcing box consis- 
tency that we just sketched can be implemented by enforcing hull consistency 
on a decomposition of the original constraint. So the procedure for enforcing 
box consistency can be seen as consisting of a number of procedures including 
the one that enforces hull consistency. One could apply the same technique to 
the arithmetic constraints on integer intervals, replacing the enforcement of hull 
consistency by one of our approaches to constraint propagation. This would 
lead to an integer equivalent of the box consistency. The efficiency of the result- 
ing procedure depends on the choice of the underlying approach to constraint 
propagation, which provides another argument for the efficiency analysis of the 
approaches here considered. 

9 A Characterization of the MULTIPLICATION 
Rules 

It is useful to reflect on the effect of the proof rules used to achieve constraint 
propagation. In this section, by way of example, we focus on the MULTIPLICA- 
TION rules and characterize their effect using the notion of bounds consistency 
as defined in [15], limited to integer intervals. Let us recall first the definition 
that we adopt here to the multiplication constraint. Given an integer interval 
[l..h] we denote by [I, h] the corresponding real interval. 

Definition 9.1 The CSP (x ■ y — z ; x G [l x ..h x ],y G [l y ..h y ],z G [l z ..h z ]} is 
called bounds consistent if 

• Va G {l x , h x } 3b G [l y , h y ] 3c G [l z ,h z ] a ■ b = c, 

• V6 G {l y , h y } 3a e [l x ,h x ] 3c G [l z ,h z ] a-b = c, 

• Vc G {l z , h z } 3a G [l x , h x ] 3b G [ly, h y ] a ■ b = a □ 

The following result entails that the MULTIPLICATION rules will not re- 
duce a CSP beyond bounds consistency. 

Theorem 9.2 (Bounds consistency) Suppose a CSP (x-y = z ; x G D x , y G 
D y ,z G D z ) with the integer interval domains is bounds consistent. Then it is 
closed under the applications of the MULTIPLICATION 1,2 and 3 rules. 
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Proof. See the Appendix. 



□ 



This theorem shows that the MULTIPLICATION rules entail a notion of 
local consistency, say M-consistency, that is implied by bounds consistency. 
However, M-consistency does not imply bounds consistency. Here is an example. 
Consider the CSP 



It is not bounds consistent, since for y — —3 no real values a G [—2, 1] and 
c G [8, 10] exist such that a ■ (—3) = c. Indeed, it is easy to check that 

{yeTZ\3xe [-2, 1] 3z e [8, 10] x ■ y = z} = (-oo, -4] U [8, oo). 

However, this CSP is closed under the applications of the MULTIPLICA- 
TION 1, 2 and 3 rules since 

• [8. .10] C int([-2..1] ■ [-3.. 10]), as int([-2..1] • [-3. .10]) = [-20.. 10], 

• [-2..1] C int([8..10]/[-3..10]) as int([8..10]/[-3..10]) = [-10. .10], and 

• [-3.. 10] C int([8..10]/[-2..1]) as int([8..10]/[-2..1]) = [-10. .10]. 

The following result clarifies that this example identifies the only cause of 
discrepancy between M-consistency and bounds consistency. Here, given an 
integer interval D := [l..h] we define (D) := {x G Z \ I < x < h}. 

Theorem 9.3 (Bounds consistency 2) Consider a CSP (f> := (x-y = z ; x G 
D x ,y G D y , z £ D z ) with non-empty integer interval domains and such that 



Suppose (f> is closed under the applications of the MULTIPLICATION 1,2 and 
3 rules. Then it is bounds consistent. 



Consequently the MULTIPLICATION rules only fail to enforce bounds con- 
sistency for the constraint x ■ y = z in case the domains of x and y are both 
of the form [l..h], with / < and h > while z can assume either only posi- 
tive numbers, or only negative numbers. Because the zeroes in the domains of 
x and y do not contribute to any solution, we can remedy this effect by tem- 
porarily splitting these domains in a positive interval and a negative interval. 
Bounds consistency for the constraint x ■ y — z is then achieved by applying 
the MULTIPLICATION rules to the resulting subproblems, and updating the 
domain of each variable with the int(-) closure of the union of its domain in 
these subproblems. 

In [2D] similar rules to our MULTIPLICATION rules are defined that apply 
this technique directly. They were defined without the use of interval arith- 
metic. It is also shown there that the LINEAR EQUALITY and LINEAR 
INEQUALITY rules enforce bounds consistency. 



(x-y = z\ xe [-2..1],ye [-3..10],ze [8.. 10]). 



S (D x ) n (D y ) implies G D z . 



(5) 



Proof. See the Appendix. 



□ 
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10 Implementation Details 



10.1 Weak Division 

We already mentioned in Section [3] that the division operation on the inter- 
vals does not admit an efficient implementation. The reason is that the int(-) 
closure of the interval division [a..6]/[c..d] requires an auxiliary computation 
in case when g" [c.d]. The preprocessing of [c.d] becomes impractical for 
small intervals [a., b], and large [c.d], occurring for example for the constraint 
n"=i x i = n"=i of the benchmark problem mentioned in Subsection ll.il This 
can be remedied by using the following variant of the division operation. We 
call it weak division since it yields a larger set (and so is 'weaker'). 

( [\min(A)] .. [max(A)\] if g [c.d], or 
[a..b] : [c.d] := < £ [a..b] and G {c, d} and c < d, 

[ [a..b]/[c..d] otherwise 

where A = {a/c',a/d',b/c',b/d'}, and [c'..d'} = [c.d] - {0}. 

Then int([a..6] : [c.d]) can be computed by a straightforward case analysis 
already used for int([a..6]/[c..d]) but now without any auxiliary computation. 
The weak division operator gives rise to the following versions of the MULTI- 
PLICATION rules 2 and 3: 

MULTIPLICATION 2w 

(x -y = z ; x £ D x , y £ D y , z £ D z ) 
(x-y = z; x G D' x ,y G D y ,z £ D z ) 

MULTIPLICATION 3w 

(x-y = z ; x £ D x , y £ D y , z £ D z ) 
(x-y = z; x £ D x , y £ D' y , z £ D z ) 

where D' x := D x n mt{D z : D y ), and D' y := D y n int(D z : D x ). 

In the assumed framework based on constraint propagation and tree search, 
all domains become eventually singletons or empty sets. It can easily be verified 
that both division operations are then equal, i.e., [a..b] : [c.d] = [a.. b] /[c.d], for 
a > b and c > d. For this reason, we can safely replace any of the reduction rules 
introduced in this paper, notably POLYNOMIAL EQUALITY, POLYNOMIAL 
INEQUALITY, and MULTIPLICATION 2 and 3, by their counterparts based 
on the weak division. For the MULTIPLICATION rules specifically, the follow- 
ing theorem states that both sets of rules actually achieve the same constraint 
propagation. 

Theorem 10.1 (MULTIPLICATION) A CSP (x ■ y = z ; x £ D x , y £ 

D y ,z £ D z ) with the integer interval domains is closed under the applications 
of the MULTIPLICATION 1, 2 and 3 rules iff it is closed under the applications 
of the MULTIPLICATION 1, 2w and 3w rules. 
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Proof. See the Appendix. 
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Let us clarify now the relation between the MULTIPLICATION rules and 
the corresponding rules based on real interval arithmetic combined with the 
rounding of the resulting real intervals inwards to the largest integer intervals. 
The CSP (x-y = z ; x G [-3..3], y G [-l..l],z G [1..2]), which we already 
discussed in the introduction, shows that these approaches yield different results. 
Indeed, using the MULTIPLICATION rule 2 we can reduce the domain of x to 
[—2.. 2], while the latter approach yields no reduction. 

On the other hand, the applications of the MULTIPLICATION rules 2w and 
3w to (x -y — z ; x £ D x ,y G D y , z G D z ) such that int(D z : D x ) ^ mt(D z /D x ) 
and int(D z : D y ) ^ int(D z / D y ) (so in cases when the use of the weak interval 
division differs from the use of the interval division) do coincide with the just 
discussed approach based on real interval arithmetic and inward rounding. This 
is a consequence of the way the multiplication and division of real intervals are 
defined, see |15j . However, as we already stated in the introduction, we believe 
that the limited precision of floating-point interval arithmetic, and the modest 
overhead of arbitrary length integers justify a separate implementation of these 
rules for arithmetic constraints on integer intervals. 

10.2 Implementation 
Platform 

Our experiments were performed using OpenSolver |23j . an experimental con- 
straint solver based on constraint propagation and tree search. OpenSolver can 
be configured by software plug-ins in a number of predefined categories, corre- 
sponding to different aspects of constraint propagation and tree search, which 
makes it particularly well-suited for carrying out comparative studies of imple- 
mentations of constraint solvers. The categories of plug-ins that are relevant for 
the experiments reported here are: 

• variable domain types, which implement the domains of variables, 

• domain reduction functions (DRFs), which correspond to the reduction 
rules, 

• schedulers of DRFs, which determine the order in which the DRFs are 
applied, 

• branching strategies, which split the search tree after constraint propaga- 
tion has terminated, and 

• several categories corresponding to different aspects of a search strategy 
that determine how to traverse a search tree. 

All experiments were performed using the Integerlnterval variable domain 
type plug-in. Domains of this type consist of an indication of the type of the 
interval (bounded, unbounded, left/right-bounded, or empty), and a pair of 
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arbitrary length integer bounds. This plug-in, and the DRFs operating on it are 
built using the already mentioned GNU MP library, which provides arbitrary 
length integers and arithmetic operations on them, including operations for 
rounding the outcome of divisions and root extractions in the desired direction. 

The branching strategy that we used selects variables using a chronological 
ordering in which the auxiliary variables come last. The domain of the selected 
variable is split into two subdomains using bisection, so the resulting search 
trees are binary trees. In all experiments we searched for all solutions, travers- 
ing the entire search tree by means of depth-first leftmost-first chronological 
backtracking. 

For the experiments in this paper a DRF plug-in has been developed that 
implements the domain reduction rules discussed in the previous sections. Every 
constraint of a CSP is enforced by a number of instantiations of this DRF: one 
for each variable occurrence. 

The scheduler plug-in that we used in the experiments maintains a flag per 
DRF, indicating whether the DRF is pending application or not. Initially, all 
DRFs are pending application. If the application of a DRF (or the branching 
strategy) modifies the domains of one or more variables, all DRFs whose output 
depends on these variables become pending application. Since in general — as 
illustrated by the example at the beginning of Section [5] — the DRFs are non- 
idempotent, this may include the DRF that has just been applied. By default, 
the scheduler plug-in keeps cycling through the set of DRFs for a given CSP 
in a specified order, applying those DRFs that are pending application. The 
cycling stops when no DRF is pending application, or when the domain of a 
variable becomes empty. 

Scheduling of Reduction Rules 

It was already shown in [22| that controlling the order in which variables are 
updated can improve the efficiency of constraint propagation algorithms, and 
for this purpose, our scheduler plug-in can be supplied with a schedule. Such 
a schedule is a sequence of indices into the set of DRFs that describes the order 
in which the scheduler will visit them, as an alternative to cycling. This is used 
in combination with full and partial decomposition, where we distinguish user 
constraints from the constraints that are introduced to define the values of 
auxiliary variables. Before considering for execution a DRF / that is part of the 
implementation of a user constraint, we make sure that all auxiliary variables 
that / relies on are updated. For this purpose, the indices of the DRFs that 
update these variables precede the index of / in the schedule. If / can change 
the value of an auxiliary variable, its index is followed by the indices of the 
DRFs that propagate back these changes to the variables that define the value 
of this auxiliary variable. 

For example, rewriting x 3 ■ y — x < 40 to simple constraints introduces an 
auxiliary variable it, which is equated with x 3 ■ y. This leads to five reduction 
rules: one for each occurrence of a variable after the rewriting step. We number 
these reduction rules as follows, where we underline in the constraint the variable 
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that is updated by the rule: 

1. u — x 3 ■ y 2. u = x 3 • y 3. u — x 3 ■ y 
4. u - x < 40 5. u - x < 40 

The fragment of the generated schedule that corresponds to enforcing the con- 
straint x 3 ■ y — x < 40 is then 4,2,3,1,5. Rules 4 and 5 correspond to the original 
inequality, but rule 4 potentially modifies u, so in the schedule, rule 4 is followed 
by rules 2 and 3, that propagate any changes to u back to x and y. Before con- 
sidering rule 5 for application, the schedule specifies that first rule 1 should be 
considered, so that any changes to the domains of x and y are propagated to 
the domain of u. 

To see that an appropriate scheduling of the rules can be beneficial compared 
to cycling through the rules, suppose that all rules are pending application, and 
that D x = D y = [1..100], and D z = Z. If we iterate the rules in their original 
order 1,2,3,4,5 then we first reduce D u to [1..100 4 ] by means of rule 1. Next, 
rules 2 and 3 are executed without making any changes. Rule 4 then reduces 
D u to [1..140], which makes rules 2 and 3 pending application again. Next, rule 
5 is executed without reducing D x . Because x and y have not changed, rule 1 is 
not set to pending application, and rule 2 is the first rule that is applied in the 
second cycle, which reduces D x to [1..5]. If use the generated schedule 4,2,3,1,5 
instead, the same reduction is achieved immediately after applying the first two 
rules, instead of the six rules that are applied if we just cycle through the rules. 

For full decomposition, there can be hierarchical dependencies between auxil- 
iary variables. Much like the HC4revise procedure of [4], the generated schedule 
then specifies a bottom-up traversal of this hierarchy in a forward evaluation 
phase, and a top-down traversal in a backward propagation phase. These phases 
are performed before and after applying a DRF of a user constraint, respectively. 
In the forward evaluation phase, the DRFs that are executed correspond to the 
MULTIPLICATION 1 and EXPONENTIATION rules. The DRFs of the back- 
ward propagation phase correspond to the MULTIPLICATION 2 and 3, and 
ROOT EXTRACTION rules. The HC4revise procedure is part of the HC4 al- 
gorithm, which enforces hull consistency for constraints on the reals using an 
implicit decomposition. For a discussion of this algorithm in the context of 
controlled constraint propagation, see [12] . 

Constraint Rewriting 

The proposed approaches were implemented by first rewriting arithmetic con- 
straints to polynomial constraints, and then to a sequence of DRFs that corre- 
spond to the rules of the approach used. We implemented the following variants: 

d u (direct, unoptimized) : the direct approach, discussed in Section where 
we isolate all variable occurrences in the original constraints without de- 
composing them first; 
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d Q (direct, optimized): the optimization of the direct approach, discussed 
at the end of Section [5l which involves dividing out common powers of 
variables in the extended arithmetic expressions that arise from isolating 
the variable occurrences; 

p u (partial, unoptimizcd): partial decomposition into simple constraints, as 
discussed in Section [SJ The decomposition is implemented by introducing 
an auxiliary variable for every nonlinear power product. This procedure 
may introduce more auxiliary variables than necessary; 

p (partial, optimized): an optimized version of variant p u , where we stop 
introducing auxiliary variables as soon as the constraints contain no more 
duplicate occurrences of variables; 

fm (full, multiplication): full decomposition into atomic constraints, as dis- 
cussed in Section allowing only linear constraints and multiplication as 
atomic constraints; 

f s (full, squaring): idem, but also allowing x = y 2 as an atomic constraint; 

f e (full, exponentiation): idem, allowing x — y n for all n > 1 as an atomic 
constraint. 

If the distinction between the different variants of an approach is irrelevant, we 
will sometimes omit the subscripts to the names d, p, and f . 

Full and partial decomposition are implemented as a rewrite step, where the 
auxiliary variables are introduced. The resulting CSP is then rewritten using 
the direct approach. During the first rewrite step the hierarchical relations 
between the auxiliary variables are recorded, and the schedules are generated 
as a part of the second rewrite step. For variants p Q and f the question of 
which auxiliary variables to introduce is an optimization problem in itself. Some 
choices result in more auxiliary variables than others. We have not treated 
this issue as an optimization problem but relied on the (somewhat arbitrary) 
heuristics described below. For this reason we have to consider the possibility 
that performance of variants p and f can be further improved because in our 
experiments we used a suboptimal decomposition. The heuristics are as follows. 

• For variant p Q we replace nonlinear power products from left to right, 
so the rightmost nonlinear term of a polynomial constraint is always left 
intact. 

• For the full decomposition approach, nonlinear power products are pro- 
cessed in the order in which they occur in the problem statement, after 
normalization to polynomial constraints. On the first occurrence of a non- 
linear power product, we start introducing auxiliary variables for terms 
that divide the power product by multiplying or exponentiating existing 
variables, and keep doing so until we have introduced an auxiliary variable 
that corresponds to the full power product. When there are several choices 
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for which existing variables to multiply or exponentiate, we introduce an 
auxiliary variable for a term with the largest possible sum of exponents, 
thereby giving preference to exponentiation over multiplication, insofar as 
it is allowed by the variant. For variant f e we first introduce auxiliary 
variables for all exponentiations in the power product. For variant f s , we 
first introduce auxiliary variables for all exponentiations that divide the 
power product, and whose exponent is a power of 2. Unused auxiliary 
variables are deleted at a later stage. 

To illustrate the latter heuristic, suppose we want to introduce an auxiliary 
variable for the term x 5 -y 3 -z. If we allow exponentiation, we start by introducing 
auxiliary variables U\ and Ui for the exponentiations in the term, and constrain 
them as follows: u\ = x 5 , 112 — y 3 ■ Next we can introduce an auxiliary variable 
U3 for x 5 ■ y 3 , x 5 ■ z, or y 3 ■ z by adding a constraint that multiplies two of the 
variables Ui, U2, and z. Because the sum of exponents is highest for the first 
option, we add U3 = u\-U2- Finally m is introduced to replace the original term: 
Ui = U3 ■ z. With only squaring allowed, we would be making these decisions 
in the presence of auxiliary variables for x 2 , a; 4 , and y 2 , where x A is obtained 
by squaring x 2 . In this case, the first auxiliary variable introduced would be 
for x 4 ■ y 2 . With only multiplication allowed, after introducing u\ = x ■ x and 
u 2 =U\'U\, we would be expanding the term be repeatedly multiplying it with 
x, y, or z. 

Except for the optimized version of the direct approach, our current imple- 
mentation can be optimized further by adopting the two-step propagation of 
linear constraints described in [14] , as discussed at the end of Section |4j Be- 
cause linear constraints are never decomposed, the effect is essentially the same 
for all alternatives that we discussed, so we have not considered this technique 
in our evaluation. 

11 Experiments 
11.1 Problems 

For evaluating the alternative approaches, we used the integer problems de- 
scribed below. Problems with only integer variables and arithmetic constraints 
are rare in practice, and in that sense, our benchmark problems are artificial, 
but they serve well to generate a purely integer workload. The approach that 
works best on these problems can also be expected to work well in a hybrid 
setting, where integer variables and arithmetic constraints are mixed with other 
types of variables and constraints. In that case, only a fraction of the workload 
will be devoted to integer arithmetic, but with the results of this study we can 
be confident that we are dealing with these constraints in an efficient way. 
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Cubes The problem is to find all natural numbers n < 10 5 that are a sum of 
four different cubes, for example 

l 3 + 2 3 + 3 3 + 4 3 = 100. 

This problem is modeled as follows: 

(1 < %1> X\ < X2 — 1, X2 < £3 — 1, X3 < £4 — 1, X4 < n, 

x\ + x\ + x 3 + x 3 = n; n S [1..10 5 ], xi,X2,x 3 ,x4 e Z) 

Opt We are interested in finding a solution to the constraint x 3 + y 2 — z 3 in 
the integer interval [L.fO 5 ] for which the value of 2x ■ y — z is maximal. 

Fractions This problem is taken from [115]: find distinct nonzero digits such 
that the following equation holds: 

A D G 
~BC + ~EF + ~HI ~ 

There is a variable for each letter. The initial domains are [1..9]. To avoid 
symmetric solutions an ordering is imposed: 

A D G 
~BC ~ ~~EF ~ ~HI 

Also two redundant constraints are added: 

A G 

3 > f and 3 < 1 

BC ~ HI ~ 

Because division is not present in our arithmetic expressions, the above con- 
straints are multiplied by the denominators of the fractions to obtain arithmetic 
constraints. We studied a representation for this problem using one equality and 
four inequalities for the ordering and the redundant constraints, and 36 dise- 
qualities A j= B, A ^ C, H ^ I. 

Kyoto The problem (see [TU]) is to find the number n such that the alphanu- 
meric equation 

KYOTO 
KYOTO 
+ K Y O T O 
TOKYO 

has a solution in the base-n number system. Our representation uses a variable 
for each letter and one variable for the base number. The variables K and T may 
not be zero. There is one large constraint for the addition, 6 disequalities K ^ Y 
... T ^ O and four constraints stating that the individual digits K, Y, O, T, are 
smaller than the base number. To spend some CPU time, we searched base 
numbers 2.. 100. 
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Sumprod This is the problem cited in Subsection ll.il for n = 14. We use the 
following representation: 

(xi + . . . + X n = Ci + . . . + Cn, 

X\ . . . x n — C\ • . . . c n , 

X\ < x 2 ,x 2 <x 3 ,.. .,x n -i < x n ; 

Xi,...,x n e [1-n], 

ci e {l},c 2 e {2},...,c„ e {n}) 

For n = 14, the value of the expression Yl7=i * ecma l s 14!, which exceeds 2 32 , 
and to avoid problems with the input of large numbers, we used bound variables 
Ci, . . . , c„ and constraint propagation to evaluate it. 

11.2 Results 

Tables [5] and [3] compare the implemented variants of our approaches on the 
problems defined in the previous subsection. The first two columns of table [2] 
list the number of variables and DRFs that were used. Column nodes lists the 
size of the search tree, including failures and solutions. The next two columns 
list the number of times that a DRF was applied, and the percentage of these 
applications that the domain of a variable was actually modified. For the opt 
problem, the DRF that implements the optimization is not counted, and its 
application is not taken into account. The reported CPU times are user time in 
seconds, as reported by the UNIX time command on a 1200 MHz Athlon CPU. 
The last column compares the performance of our implementation to that of 
ECL 4 PS e , and will be discussed at the end of this section. 

Table [3] lists measured numbers of basic interval operations. Note that for 
variant d Q , there are two versions of the division and addition operations: one 
for integer intervals, and one for intervals of reals of which the bounds are 
rational numbers (marked Q). Columns multl and multF list the numbers of 
multiplications of two integer intervals, and of an integer interval and an integer 
factor, respectively. These are different operations in our implementation. 

For the cubes, opt, and sumprod problems, the constraints are already in 
simple form, so variants d u , d Q and p are identical. For cubes and opt all 
nonlinear terms involve a single multiplication or exponentiation, so for these 
experiments also variants p u and f e are the same. For the fractions problem, 
and for sumprod, no exponentiations are used, so all three variants of the full 
decomposition approach that we implemented are identical. 

The results of these experiments clearly show the disadvantage of implement- 
ing exponentiation by means of multiplication: there is less domain reduction 
because we increase the number of variable occurrences (see the dependency 
problem, discussed in Section [8]). For opt and variant f m , the run did not 
complete within reasonable time and was aborted. 

For fractions the symbolic manipulation of variant d Q reduces the search tree 
by a factor 0.70. However, this reduction is not reflected in the timings, and in 
fact the CPU time even increases. The reason is that computing the domain 
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DRFs 




CPU 






nvar 


nDRF 


nodes 


applied 


%eff. 


(sec.) 


EC17PS e 


cubes 














6.54s 


d,p 


5 


14 


169,755 


1,876,192 


9.52 


9.69 




Pu,fe 


9 


22 


169,755 


2,237,590 


16.28 


6.53 





fm 


13 


34 


206,405 


3,011,749 


20.02 


8.53 




fs 


13 


34 


178,781 


2,895,717 


20.62 


8.80 




opt 














5752.70s 


d,p 


4 


7 


115,469 


5,187,002 


42.16 


21.55 


+ 


Pu,fe 


8 


15 


115,469 


9,800,017 


60.00 


22.75 


+ 


fm 


10 


21 


? 


? 


7 


7 




fs 


10 


21 


5,065,195 


156,906,444 


46.49 


422.93 




fractions 














6.90s 


d u 


9 


154 


11,289 


1,193,579 


3.65 


15.40 




d 


9 


154 


7,879 


734,980 


3.45 


17.38 




Pu 


37 


210 


11,289 


1,410,436 


23.27 


4.89 




Po 


32 


200 


11,289 


1,385,933 


21.65 


5.25 




f 


43 


208 


11,131 


1,426,204 


27.76 


4.98 




kyoto 
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d u 


5 


37 


87,085 


3,299,814 


6.09 


21.84 




do 


5 


37 


87,085 


3,288,461 


5.94 


44.56 


+ 


Pu 


13 
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87,085 


3,781,514 


23.02 


10.93 




Po 


12 


51 


87,085 


3,622,461 


21.45 


11.24 




fm 
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60 


87,087 


4,276,066 


26.70 


10.40 




fs 


16 
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87,085 


4,275,957 
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10.39 




fe 


16 


59 


87,085 


3,746,532 


23.26 


9.42 
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23.25s 


d,p 


28 


82 


230,233 


10,910,441 


7.91 


102.49 




Pu 


30 


86 


230,233 


9,196,772 


9.39 


80.59 




f 


54 


134 


55,385 


3,078,649 


18.01 


23.75 





Table 2: Statistics and comparison with ECL J PS e 
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4,245 


14,408 
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f S 
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4,305 
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2,299 


4,599 


1,443 


1,444 
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26,037 


Pu,fe 


1,636 


1,538 


2,150 


738 


8,138 


4,445 
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fm 
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? 
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? 


? 
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21,066 


18,106 
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29,584 








1,550 Q 






1,355 Q 
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734 


933 


4,736 


4,669 


11,071 


Po 








776 


1,509 


5,292 


5,147 


12,725 


f 








693 


339 


4,835 


4,769 


10,636 
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d u 


735 


11,041 


1,963 


13,853 


10,853 


13,946 


52,390 


d 


735 


8,146 


218 


8,955 


12,516 


10,592 


48,749 








4,310 Q 






3,277 Q 




Pu 


383 


759 


1,591 


484 


5,324 


7,504 


16,044 


Po 


383 


759 


1,597 


1,360 


5,756 


8,008 


17,863 


fm 








1,991 


578 


5,324 


7,505 


15,398 


fs 


< 0.5 


< 0.5 


1,990 


578 


5,324 


7,504 


15,397 


fe 


1 


1 


1,554 


484 


5,324 


7,504 


14,868 
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d,p 








4,032 


100,791 


85,419 


149,479 


339,721 


Pu 








2,186 


27,948 


81,728 


149,479 


261,340 


f 








609 


205 


25,799 


46,960 


73,573 



Table 3: Measured numbers (thousands) of interval operations 
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updates involves adding intervals of real numbers. The arithmetic operations on 
such intervals are more expensive than their counterparts on integer intervals, 
because the bounds have to be maintained as rational numbers. Arithmetic 
operations on rational numbers are more expensive because they involve the 
computation of greatest common divisors. For kyoto the symbolic manipulation 
did not reduce the size of the search tree, so the effect is even more severe. 

In general, the introduction of auxiliary variables leads to a reduction of the 
number of interval operations compared to the direct approach. As discussed at 
the end of Section [51 this is because auxiliary variables prevent the evaluation 
of subexpressions that did not change. This effect is strongest for fractions, 
where the main constraint contains a large number of different power prod- 
ucts. Without auxiliary variables all power products are evaluated for every 
POLYNOMIAL EQUALITY rule defined by this constraint, even those power 
products the variable domains of which did not change. With auxiliary vari- 
ables the intervals for such unmodified terms are available immediately, which 
leads to a significant reduction of the number of interval multiplications. For 
sumprod, the difference between variants d and p u is a bit artificial, because 
the operations that are saved involve the computation of the constant term 
C\ ■ . . . ■ c n . A comparable number of interval additions can be saved if we in- 
troduce a variable for the constant term C\ + . . . + c n . If we add these variables 
to the CSP all variants of the direct and partial decomposition approaches are 
essentially the same. 

That stronger reduction is achieved as a result of full decomposition, men- 
tioned in Section [3 is seen for the fractions benchmark and more prominently 
for sumprod. In the latter benchmark, this effect depends on a decomposition 
of the term n"=i 

X% clS X 1 • (x2 •(•■■• {x n -i ■ x n ) . . .)), with an auxiliary vari- 
able for each pair of matching brackets. The decomposition then matches the 
chronological ordering used to select the variable for branching. If the ordering 
is reversed, the number of nodes is equal to that of the other approaches. The 
effect described in Section [5] is not demonstrated by these experiments. 

If we do not consider the symbolic manipulation of variant d Q , variant f e 
leads to the smallest total number of interval operations in all cases, but the 
scheduling mechanism discussed in Section [TU] is essential for a consistent good 
performance. If for example the schedule is omitted for opt, the number of 
interval operations almost triples, and performance of variants p u and f e is then 
much worse than that of d u . This is conform the observations of [12], where 
it is demonstrated that for constraints on reals, enforcing hull consistency for a 
decomposition through repeated application of the HC4revise procedure yields 
superior performance compared to the basic HC3 algorithm. Based on these 
observations, we expect that the benefit of using the schedule will grow with 
the number of variables. 

The total numbers of interval operations in table [3] do not fully explain all 
differences in elapsed times. One of the reasons is that different interval opera- 
tions have different costs. Also some overhead is involved in applying a reduction 
rule, so if the number of applications differs significantly for two experiments, 
this influences the elapsed times as well {opt, d, p u ). The elapsed times are 
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not the only measure that is subject to implementation details. For example, 
we implemented division by a constant interval [— 1.. — 1] as multiplication by 
a constant, which is more efficient in our implementation. Such decisions are 
reflected in the numbers reported in table [31 

For each of the benchmarks, the last column of Table [5] compares the per- 
formance of the variants that we implemented with that of an ECL J PS e [7] 
program that directly encodes the problem statement of Subsection 111.11 using 
the ic library. For each problem, the first entry in this column lists the CPU 
time reported by ECL'PS 6 for an all-solution search, where we applied the same 
branching scheme as we used in OpenSolver. The other entries compare prop- 
agation strength, for which we ran the solvers without search, and compared 
the resulting domains of the variables. A mark '=' means that the computed 
domains are the same, '+' that our variant achieved stronger reduction, and '-' 
that constraint propagation is weaker than with ECL/PS 6 . 

In addition, for cubes we verified that the number of nodes in the ECL 4 PS e 
search tree is identical to that for all variants except f m and f s , which nicely 
fits with the comparable CPU times. In contrast, for the kyoto benchmark, the 
number of nodes in the search tree is substantially lower for our approaches than 
for ECI/PS e , and so is the CPU time. For the opt problem the CPU time for 
our approaches (except f m ) is also substantially lower than for ECL l PS e . We 
have not verified the number of nodes visited by the minimize/2 built-in, but 
the sequence of suboptimal solutions is identical to that found by our approaches 
(not verified for f m ). For this comparison we used ECL'PS 6 version 5.10. 

12 Conclusions 

In this paper we discussed a number of approaches to constraint propagation 
for arithmetic constraints on integer intervals. To assess them we implemented 
them using the OpenSolver framework of [23) , and compared their performance 
on a number of benchmark problems. We can conclude that: 

• Implementation of exponentiation by multiplication gives weak reduction. 
In the full decomposition approach x = y n should be used as an atomic 
constraint. 

• The optimization of the direct approach, where common powers of vari- 
ables are divided out, can significantly reduce the size of the search tree, 
but the resulting reduction steps rely heavily on the division and addi- 
tion of rational numbers. These operations are more expensive than their 
integer counterparts, because they involve the computation of greatest 
common divisors. As a result, our implementation of this approach was 
inefficient. 

• Introducing auxiliary variables can be beneficial in two ways: it may 
strengthen constraint propagation, as discussed in Sections [6] and [3 and 
it may prevent the evaluation of subexpressions the variable domains of 
which did not change. 
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• As a result, given an appropriate scheduling of the rules, the full and par- 
tial decomposition approaches perform better than the direct approach 
without the optimization, in terms of numbers of interval operations. Ac- 
tual performance depends on many implementation aspects. However for 
our test problems the performance of variants p u , p and f e does not differ 
much, except for one case where the decomposition of a single multiplica- 
tion of all variables significantly reduced the size of the search tree. 

Because of the inherent simplicity of the reduction rules and the potential 
reduction of the search tree, full decomposition of arithmetic constraints into 
multiplication, exponentiation, and linear constraints is our method of choice. 
However, a hierarchical scheduling of the resulting reduction rules is essential 
for efficient constraint propagation, and if a solver does not provide facilities for 
controlling the propagation order, the direct approach is preferable. 

Given that the optimization of the direct approach can achieve a significant 
reduction of the search tree, it would be interesting to combine it with full de- 
composition. Depending on the effect of the symbolic manipulation, a selection 
of the optimized rules that enforce a particular constraint according to variant 
d Q could be used as redundant rules. In this case, the internal computations 
need not be precise, and we could maintain the rational bounds as floating-point 
numbers, thus avoiding the expensive computation of greatest common divisors. 
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Appendix 

We provide here the proofs of the Bounds consistency Theorems 19.21 and 19.31 
and the MULTIPLICATION Theorem [TOTj 

Proof of the Bounds consistency Theorem 19.21 

Let <j) := (x ■ y = z ; x £ D X ,D £ D y ,z £ D z ). Call a variable u of <f> 
bounds consistent if the bounds of its domain satisfy the condition of the bounds 
consistency (see Definition 19. lj) . 

Given an integer interval [l..h] denote by [l..h] the corresponding real interval 
[l,h]. Suppose that D x = [l x ..h x ], D y = [l y ..h y ], D z = [l z ..h z ]. To show that <f> 
is closed under the applications of the MULTIPLICATION 1 rule it suffices to 
prove that 

{l z ,h z }Cmt(D x -D y ). (6) 

So take c £ {l z ,h z }. By the bounds consistency of z we have c = a ■ b 
for some a £ D x and b £ D y . Since D x and D y are integer intervals we have 
[a\, \a] £ D x and \b\, [6] £ D y . To prove ([6|), by the definition of D x ■ D y , we 
need to find a\ , ai £ D x and b\ , bi £ D y such that 

a i • b\ < c < a 2 • b 2 - 

The choice of a\, a 2 , b\ and b 2 depends on the sign of a and of b and is provided 
in the following table: 



condition 


a 1 


bi 


ai 


b 2 


a = 


a 


[b\ 


a 


Vb\ 


6 = 


[a\ 


b 


[a\ 


6 


a > 0,6 > 


[a\ 


[b\ 


\a] 




a > 0,6 < 


\a] 


[b\ 


[a\ 




a < 0,6 > 


[a\ 


\b] 


\a] 


\P\ 


a < 0,6 < 


\a] 


\b] 


[a\ 


Vb\ 
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To prove that </> is closed under the applications of the MULTIPLICATION 
2 and 3 rules it suffices to prove 

{l x , h x } C mt{D z /D y ) and {l y , h y } C mt{D z /D x ). (7) 

We need to distinguish a number of cases. The case analysis depends on the 
position of w.r.t. each of the intervals D x and D y . This leads to 9 cases, which 
by symmetry between x and y can be reduced to 6 cases. We present here the 
proofs for representative 3 cases. 

Case 1. l x > 0, l y > 0. 

By the bounds consistency of x for some b £ [l y , h y ] we have l x ■ b £ [l z , h z ]. 
Then b < h y and l x > 0, so l x ■ b < l x ■ h y . Also l z < l x ■ 6, so 

h < lx ■ h y . 

Next, by the bounds consistency of y for some a £ [l x , h x ] we have a ■ h y e 
[l z ,h z ]. Then l x < a and h y > 0, so l x ■ h y < a ■ h y . Also a ■ h y < h z , so 

lx ' h y ^ h z . 

So l x ■ h y € [Z z --/iz] an d consequently by the definition of the integer intervals 
division 

lx € D z /D y and /i y e D z /D x . 
By a symmetric argument 

€ D z /D y and Z y € D z /D x . 



Case 2. l x > 0, /i y < 0. 

By the bounds consistency of x for some b e we have h x ■ b e [Z z , /i z ]. 

Then b < h y and h x > 0, so h x ■ b < h x ■ h y . Also l z <h x ■ b, so 

Next, by the bounds consistency of y for some a £ [l x , h x ] we have a ■ h y £ 
[l z ,h z ]. Then a < h x and h y < 0, so a ■ h y > h x ■ h y . Also h z > a ■ h y , so 

h x - h y < h z . 

So h x ■ h y £ [l z ..h z ] and consequently by the definition of the integer intervals 
division 

frx G D z /D y and /i y e D z /D x . 

Further, by the bounds consistency of x for some b £ [l y , h y ] we have l x ■ b £ 
[l z ,h z ]. Then l y < b and l x > 0, so l x ■ l y < l x ■ b. Also l x ■ b < h z , so 

lx ' ly — h z . 
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Next, by the bounds consistency of y for some n£ [l x , h x ] we have a-l y G [l z , h z ]. 
Then l x < a and l y < 0, so l x ■ l y > a ■ l y . Also a ■ l y > l z , so 

l Z lx ' ly- 

So l x ■ l y G [l z ..h z ] and consequently by the definition of the integer intervals 
division 

k G £>z/-Di/ and l y G £> Z /-D X . 

Case 5. l x <0 < h x , l y >0. 

The proof for this case is somewhat more elaborate. By the bounds consis- 
tency of x for some b G [/ y , /i y ] we have l x ■ b £ [l z , h z \. Then l y < b and ^ < 0, 
so l x ■ l y > l x ■ b. But also l x ■ b> l z , so 

/z ^ /x ' /y ■ 

Next, by the bounds consistency of y for some a G [l x , h x ] we have a ■ l y G 
[l z ,h z ]. Then l x < a and l y > 0, so l x ■ l y < a ■ l y . But also a ■ l y < h z , so 

lx ' ly — h z • 

So l x ■ l y £ [l z ..h z ] and consequently by the definition of the integer intervals 
division 

k G Dz/Dj, and l y G 5 2 /-Dx. 

Further, by the bounds consistency of a; for some 6 G [/ y , /i y ] we have h x • b G 
[l z ,h z ]. Then l y < b and /i x > 0, so h x ■ l y < h x ■ b. But also h x ■ b < h z , so 

h x -l y < h z . 

Next, we already noted that by the bounds consistency of y for some a G 
[Z x , /i x ] we have a ■ l y G [Z z , /i z ]. Then a < h x and / y > 0, so ct • ly ^ h x • l y . But 
also l z < a- ly, so 

/ Z /l^; ' /y. 

So /la; • ly G [Z z ../i z ] and consequently by the definition of the integer intervals 
division 

h x G D z /D y . 

It remains to prove that h y G D z /D x . We showed already l x ■ l y < h z . 
Moreover, l x < and l y < h y , so l x ■ h y < l x ■ l y and hence 

^a; ' /^y ^ /^z • 

Also we showed already l z < h x ■ l y . Moreover h x > and Z y < /i y , so 
h x ■ l y < h x ■ h y and hence 

lz ^ h x • hy. 

So if either l z < l x ■ h y or h x ■ h y < h z , then either l x ■ h y G [l z ..h z ] or h x ■ h y £ 
[l z ..h z ] and consequently /i y G D z /D x . 
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If both l x ■ h y < l z and h z < h x ■ h y , then 

[lz-h z ] C [l x ..h x ] ■ hy. 

In particular for some a 6 D x we have l z = a ■ hy, so ft y € D z /D Xl as well. 
This concludes the proof for this case. □ 

Proof of the Bounds consistency Theorem 19.31 

We consider each variable in turn. We begin with x. Suppose that D x = [l x ..h x \. 
cj) is closed under the applications of the MULTIPLICATION 2 rule, so 

{l x ,h x }Cmt{D z /D y ). (8) 

To show the bounds consistency of x amounts to showing 

{lx,h x }<ZLT z (Z)LT y . (9) 

(Recall that given real intervals X and Y we denote by X Y their division, 
defined in Section [51) 

Case 1. int(D z /D y ) = Z. 

This implies that £ D z n D y , so by the definition of real intervals division 
TT z (Z)LTy = (-00,00). Hence © holds. 

Case 2. mt(D z /D y ) ^ Z. 

So mt{D z / D v ) is an integer interval, say mt(D z /D y ) = [l zy ..h zy ]. Two 
subcases arise. 

Subcase 1. D z D y is a, possibly open ended, real interval. 
By (jHJ for some 61 , 62 £ D y and ci , C2 € -D z we have 

hy ■ b\ = Ci, 

h zy ■ b 2 = c 2 . 

Let 

b := min(&i, 62), b := max(6i, 6 2 ), c := min(ci, c 2 ), c := max(ci, C2). 

So {l zy , h zy } C [c, c]0[6, 6]. Also [c, c]0[6, 6] C LT z cSlTy. Hence {Z zy , ft^} CZ^0 

D y and consequently, by the assumption for this subcase, [l zy , h zy ] C D z D^. 
This proves ([9]) since by |(5J) C [l zy ,h zy ]. 

Subcase 2. D z D y is not a, possibly open ended, real interval. 

In what follows for an integer interval D := [l..h] we write D > if I > 0, 
D < if h< 0. Also recall that (D) := {1 e Z | ( < K ft}. 

This subcase can arise only when D z > and £ (-D y ) or £) 2 < and 
G (-Dj,), see [18] (reported as Theorem 4.8 in p~5]), where the definition of the 
division of real intervals is considered. 

Since (j> is closed under the MULTIPLICATION rule 3 

D y C intp./A,). 
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So \at{D z /D x ) ^ since by assumption D y is non-empty. Also, since G" D z , 
we have mk(D z / D x ) ^ Z. So int(D z / D x ) is a non-empty integer interval such 
that G (mt(D z /D x )}. 

But D z > or D z < 0, so if D x > 0, then int(D z /D x ) > or mt(I3 z /Z> x ) < 
and if D x < 0, then int(D z /D x ) > or mt(D z /D x ) < 0, as well. So G (D x ). 
Hence G (D x ) n (A,) while ^ D 2 . This contradicts ©. So this subcase 
cannot arise. 

The proof for the variable y is symmetric to the one for the variable x. 

Consider now the variable z. </> is closed under the applications of the MUL- 
TIPLICATION 1 rule, so 

D z C mt{D x ■ D y ). 

Take now c G D z . Then there exist a\,a 2 G D x and 61,62 G -D a such that 
ai • 61 < c < <i2 • 62. We can assume that both inequalities are strict, that is, 

ai ■ bi < c < a 2 ■ 6 2 , (10) 

since otherwise the desired conclusion is established. 
Let 

a := min(ai, a 2 ), a := max(ai, a 2 ), b := min(foi, b 2 ) 1 6 :— max(&i, b 2 )- 

We now show that a G [a.. a) and b G [6. .6] exist such that c = a • b. Since 
[a.. a] C Dj, and [6. .6] C D y , this will establish the bounds consistency of z. 

The choice of a and 6 depends on the signs of a\ and 62. When one of these 
values is zero, the choice is provided in the following table, where in each case 
on the account of (fTU|) no division by zero takes place: 



condition 


a 


b 


ai = 


c/6 2 


62 


a 2 = 


c/h 


61 


61 = 


a 2 


c/a 2 


6 2 = 


ai 


c/ai 



It is straightforward to show that in each case the quotient belongs to the 
corresponding interval. For example, when a\ — we need to prove that c/62 G 
[a..a]. By (JTUJ) a 2 ^ 0. If a 2 > 0, then again by |10|), 6 2 > 0, so c/b 2 G [0..a 2 ]. 
In turn, if a 2 < 0, then also by (fT0|) 62 < 0, so, yet again by (fT0|) . c/62 G [a2..0]. 

When neither ai nor b 2 is zero, the choice of a and 6 has to be argued case 
by case. 

Case 1. ai > 0, 6 2 > 0. 

Then by (fT0|) 61 < c/a\ and c/62 < a 2 . Suppose that both b 2 < cja\ and 
c/62 < Oi< Then ai • 62 < c < ai • b 2 , which is a contradiction. So either 
c/ai < b 2 or a\ < c/b 2 , that is either c/a\ G [61. .b 2 ] or c/62 G [ai..oa]. 
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Case 2. a x > 0, 6 2 < 0. 

Then by (fT0|) b\ < c/a\ and a 2 < c/b 2 . Suppose that both 62 < c/a\ and 
ai < c/62. Then a\ ■ b 2 < c < a± ■ b 2 , which is a contradiction. So either 
c/a% < b 2 or c/b 2 < a 2 , that is either c/a\ £ [61. .62] or c/62 € [fl2..ai]. 

Case 3. a ± < 0, b 2 > 0. 

Then by (fT0|) c/ai < 61 and c/62 < a 2 . Suppose that both c/a\ < b 2 and 
c/62 < Q>i- Then a\ ■ b 2 < c < a\ ■ b 2 , which is a contradiction. So either 
62 < c/ai or a\ < c/b 2 , that is either c/a\ £ [62. .61] or c/62 £ [a\..a 2 ]. 

Case 4- a i < 0, 62 < 0. 

Then by (fTP)) c/a\ < 61 and a 2 < c/b 2 . Suppose that both c/a\ < b 2 and 
ai < c/62. Then a\ ■ b 2 < c < a\ ■ b 2 , which is a contradiction. So either 
62 < c/a% or c/62 < ai, that is either cja\ £ [62. .61] or c/62 € [02. .ai]. 

So in each of the four cases we can choose either a := a\ and 6 :— c/a\ or 
a := c/b 2 and 6 :— b 2 . □ 

Proof of the MULTIPLICATION Theorem \10A\ 

The weak interval division produces larger sets than the interval division. As a 
result the MULTIPLICATION rules 2w and 3w yield a weaker reduction than 
the original MULTIPLICATION rules 2 and 3. So it suffices to prove that 
<f> := (x ■ y = z ; x £ D x ,y £ D yi z £ D z ) is closed under the applications 
of the MULTIPLICATION 1, 2 and 3 rules assuming that it is closed under 
the applications of the MULTIPLICATION 1, 2w and 3w rules. Suppose that 
D x = [l x ..h x ],D y — [ly..h y ] 1 D z — [I g . .h g] . The assumption implies 

{l x , h x } C mt(D z : D y ) (11) 

and 

{ly,h y } Cint(A, :D X ) (12) 

The proof is by contradiction. Assume that (fTTj) and (|12p hold, while (f> 
is not closed under application of MULTIPLICATION 2 and 3. Without loss 
of generality, suppose that MULTIPLICATION 2 is the rule that can make a 
further reduction. This is the case iff 

mt(D z /D y ) C int(D z : D v ). 

By definition, the proper inclusion implies that l v > or h y < 0. Assume 
l y > 0, the case for h y < is similar. Let l' y := max(l,^), and let A :— 
{lz/ly,lz/hy,h z /l y ,h z /h y }, and B := {l z /l x ,l z /h x ,h z /l x ,h z /h x }. A further 
implication of the proper inclusion is that one or both of l' y and h y do not have 
a multiple in D z : otherwise min(A) and max(A) would be elements of D z /D y , 
and we would have int(Z? z : D y ) — int(D z / D y ). The cases for I' and h y can be 
seen in isolation, and their proofs are similar, so here we only consider the case 
that l' y does not have a multiple in D z . In what follows we can assume ^ D z , 
since otherwise I' and h y do have a multiple in D z . 
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Case 1. l z > 0. 

From (fTTjl it follows that h x < [max(y4)J, which for the case l' y ,h y ,l z ,h z > 
that we consider here implies h x < [h z /l' y \. Because [l z ..h z ] does not contain a 
multiple oily, we have \ h z /l' y \ — [l z /l y \, so 

h x < iyi' y \. 

A further consequence of (fTT]) is that l x ,h x > 0. From (|T^)l it follows that 
l' y > r mm (^)l * which for l x , l z > implies 

l'y> \lz/k x ] >l Z /h x >l z /\h z /l' y \. 

Because V is no divisor of l z , and both numbers are positive, we have \ l z /l' y \ < 
l z /l' y , and consequently l' y > l z /(l z /l' y ), leading to l' y > l' y , which is a contradic- 
tion. 

Case 2. h z < 0. 

Similarly, because l' y ,h y > and l z ,h z < 0, it follows from pip that l x > 
\mm(A)~\ = \l z /l' y ~\, and l x , h x < 0. Because [l z --lh] does not contain a multiple 
of l' y , we have \l z /l y ~\ = \h z /l' y \, so 

l x > \h z /l' y ~}. 

We use this information in the following implication of (|12|1 : 

l' y > rmin(B)] = \h z /Q > h z /l' x 

to get l' y > h z /\h z /l' y ]. Because \ \h z /l' y ]\ < \h z /l' y \, we have l' y > h z /(h z /l' y ), 
leading to I' > I' which is a contradiction. □ 
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